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A GENERALIZED THEORY FOR THE DESIGN OF CONTRACTION CONES 
AND OTHER LOW- SPEED DUCTS 

By Raymond L. Barger and John T. Bowen 
Langley Research Center 

SUMMARY 

A generalization of the Tsien method of contraction- cone design is described. The 
design velocity distribution is expressed in such a form that the required high-order 
derivatives can be obtained by recursion rather than by numerical or analytic differen- 
tiation. The method is applicable to the design of diffusers and converging-diverging 
ducts as well as contraction cones. The computer program is described and a FORTRAN 
listing of the program is provided. 


INTRODUCTION 

For incompressible flows in ducts of slowly varying radius the one-dimensional 
flow relation between the velocity and the cross-sectional area can be used to predict the 
velocity distribution in a given duct or to design a duct for a desired velocity distribution. 

However certain applications, such as the contraction cone for a wind tunnel, require 
short ducts with a relatively rapid variation of the wall radius. For such applications the 
one-dimensional relation no longer suffices, and a solution of the differential equation of 
the flow must be sought. Tsien (ref. 1) derived a solution for the stream function in terms 
of a prescribed axial velocity distribution, and applied this solution in the design of a 
wind-tunnel contraction cone. 

It should be mentioned that such a design, obtained from the incompressible-flow 
equations, is a conservative design in the sense that when it is operated at off-design 
compressible conditions, the ratio of exit velocity to entering velocity is higher than the 
design ratio. 

The major difficulty in applying Tsien' s solution arises from the stringent require- 
ments on the form of the input axial velocity distribution. These requirements are such 
that only a small class of functions can be used to describe the design velocity distribu- 
tion. This problem has been considered in reference 2, where a form of velocity distri- 
bution different from that of reference 1 is used. 



This form allows more freedom in shaping the design velocity function, but it intro- 
duces the problem of requiring hand calculations in analytic form of many derivatives of 
the function. Other analyses have utilized different formulations of the solution of the 
differential equation according to the way the variables are separated in solving the equa- 
tion. The method of reference 3 assigns an exponential type of variation in the axial 
direction, and so is limited to a single design velocity distribution. References 4 and 5 
use a periodic axial velocity distribution, but inasmuch as the flow is not periodic, this 
formulation gives rise to errors near the beginning and the exit of the contraction cone. 

An additional problem with this latter method is that the finite-term trigonometric approx- 
imation to a function is in general a function that oscillates about the desired function, 
and such a "wavy" distribution is not very satisfactory for design purposes. The three 
forms of solution for the stream function are given explicitly in reference 5. 

The present analysis represents a generalization of the method of reference 1, in 
that a wide range of design velocity distributions is permitted so that the method is no 
longer restricted to a specific contraction cone but may be applied to the design of a wide 
variety of ducts. Greater accuracy is obtained through the use of an electronic computer 
and by retaining a large number of terms in the series solution. 

SYMBOLS 

A,B,c,d n arbitrary parameters and coefficients in expression for design velocity 
f d = f g- f P 

f general form of design velocity distribution at centerline 

& 

fp preliminary form of design velocity distribution at centerline (eq. (3)) 

fg design velocity distribution at centerline used in reference 1 

G total velocity 

H n nth Hermite polynomial 

k upper summation index 

m,n indices 

x,r cylindrical coordinates 
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u,v axial and radial velocity components, respectively 

z = cx 

6 mn Kronecker delta 

_ c _ C 2 X 2 

$ = e c x 

ft 

stream function 

Subscripts: 
i initial 

f final 

m,n indices 


ANALYSIS 


Tsien's solution (ref. 1) for the axial and radial velocities in incompressible axi- 
symmetric flow is 


u 


0 2 2 V) 2 


= 1 


/ 2n-l ( 211 - 1 ), \ 

(-1) 2nr f Q ' (x) 

2 2n ( n'-) 2 


where fg( x ) is the prescribed velocity on the axis. The stream function can be obtained 
by integration. Its k-term approximation is 



(-l) n “ 1 f(^n-2)^ r 2n 
2 2n_1 n[(n - l)’.] 2 


( 1 ) 


The kind of functions fg(x) which are appropriate for describing the axial velocity 
distribution will now be examined. It is apparent that if the series is truncated at the 
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nth term 2n - 2 derivatives of fg are required. Therefore, fg must be such that 
these derivatives can be obtained in analytic form because it is generally impossible to 
obtain high-order derivatives numerically with accuracy. Furthermore, as has been 
pointed out in reference 2, the simplest way to insure that conditions at infinity upstream 
and downstream be uniform is to require that all the derivatives of fQ vanish as x — ±oo } 
but of course fg must not itself vanish at x = ±°°. 

Thus it is seen that the class of functions that can be used to describe the axial 
velocity is severely limited. 

In reference 1 this velocity distribution is prescribed by the function 

Vo< x > = f o< x > = °- 65 + ^.C e " 2dx (2> 

The following analysis generalizes the procedure of reference 1 so that a much larger 
variety of design velocity distributions is permitted and in such a way that the series can 
be carried out to an arbitrary number of terms without any penalty except a trivial 
increase in machine computing time. 

The basic preliminary form of the design velocity function is chosen to be 


f P (x,=A + 2 f lo e ’ c2x2dx 


(3) 


which is only a slight generalization of equation (2). The derivatives of this expression 
can be obtained by recursion, following a development similar to that of reference 1: Let 


. . c -c^x^ 

<f.(x) = e c x 

V7T 

then the m + 1 derivative of f p is 
4 m+1) (x) = 2B 4> (m) (x) 


Substitute 



then 


( 4 ) 


$( m >(x) = -^c m -^-e- z = 


_m+l „ _2 

c -(- l) m e“ z H m (z) 


dz 


4 * 
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or 


$(m){x) _ (-ii) m c m ^(x) H m (z) (5) 

where H m (z) is the mth Hermite polynomial (see eq. (29) on p. 91 of ref. 6). The 
recurrence formula for Hermite polynomials is 

H m (z) = 2z H m _ 1 (z) - 2(m - 1) H m _ 2 (z) 

Multiply both sides by (-l) m c m 3>(x) and then equation (5) becomes 

* (m) (x) = 2z(-l) m c m H m _ 1 (z) - 2(m - l)(-l) m c m H m _ 2 (z) 

= -2 C 2 [(- 1) 111- 1 |-c m- * H m _ 1 (z) + (m - l)(-l) m - 2 c m - 2 H m _ 2 (z) 

= -2c2|x $( m- l)(x) + (m - 1) $ (m_2) (x)] 

which is the desired recurrence formula for the derivatives. 

Now consider a more general design velocity function f g (x), obtained by adding 
terms to f p (x). Since the initial and final velocities are determined by the coefficients 
in fp(x), these additional terms and all their derivatives must vanish at ±°o. They should 
also be such that an arbitrary number of differentiations can be performed analytically 
in a simple manner. These conditions are satisfied by the form: 

k 

f g(x) = f p( x ) + Yj d n e " x2 H nW ( 6 ) 

0 

2 

The factor e A in each term of the series assures that the conditions at ±°° will not 
be affected. The derivatives of these terms are obtained by the recurrence formula: 

^_e’ x2 H n (x)_ = -e' x2 H n+1 (x) 

(ref. 7, p. 786, where the stated formula contains an extraneous factor of 2). 

The coefficients d n can be determined as follows by means of the orthogonality 
property of the Hermite polynomials. Denoting 

f d (x) = f g (x) - f p (x) 
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and substituting in equation (6) 


IS. 

f d (x) =^d n e‘ x2 H n (x) 
0 


( 7 ) 


multiplying by H m (x) and integrating, yields 


oo k 2 r 

J fd(x) H m (x) ^ = ^ d n J e ’ X H n( x ) H m( x ) ** = J d n 

n=0 


k 

E 

n=0 


Thus 


d m = -S i 1 " f f d( x ) H m( x ) 

2 m'.^TT^-oo 

where the finiteness of the integral is assured by the nature of f d (x) . Of course the 
series in equation (7) will, in general, only approximate f d (x) since only k terms are 
used. 

A simpler, but less accurate, approximation could be obtained by matching the func- 
tion f d (x) at k points and obtaining the coefficients as solutions of k simultaneous 
linear equations. 

Actually, neither of these methods for determining the coefficients has been used 
so far. Rather, the calculation of fg(x) was programed for visual display, and various 
values of the coefficients were tried until a close approximation to the desired fg(x) was 
obtained. 

A description and listing of the computer program is given in the appendix. 


DESIGN PROCEDURE 


The basic considerations that govern the design of a contraction cone from a pre- 
scribed axial velocity distribution have been discussed in reference 2. In general the 
same considerations are applicable to the design of other kinds of ducts. 

After selecting an appropriate axial velocity function the next step in the procedure 
is to determine several streamlines by solving the equation for the stream function, 

^(-l) n_1 4 2n - 2) (x) r 2n 
1 2 2n_1 n[(n - l)’] 2 
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where 



for r at the designated x- stations with fixed values of \f/. A computer program library 
routine, utilizing interval-halving, was used for this purpose. It may be noted that, in 
accordance with Descarte's rule of signs, there may be as many as k positive solutions 
of equation (1) for r2 and so for r (for fixed and x). However, any possible 
ambiguity in the solution can be avoided by making an initial estimate of the radius from 
the one-dimensional approximate relation between velocity and area ratio, after one point 
on the streamline has been computed. 

As successive streamlines are determined, the velocity distributions along the 
streamlines are also computed. These display a greater radial variation in regions of 
larger curvature, eventually leading to an adverse velocity gradient in regions of inward 
turning of the wall. Of course some radial velocity gradient is normally acceptable, and 
generally a slight adverse velocity gradient can be tolerated by the boundary layer. These 
factors must be considered when selecting a streamline for the actual duct contour inas- 
much as the duct length is shortened by taking larger values of the stream function. Since 
a short duct implies savings in material, space, and wall-friction losses, the usual design 
goal is to have the shortest possible duct compatible with acceptable flow quality. 

DISCUSSION AND EXAMPLES 


is determined by the choice of the vari- 
can be readily demonstrated. Using the 

that upstream, at x = -°°, 

f g)i = A- B 

and downstream, at x = +«>, 
f g,f = A + B 


The form of the design velocity distribution 
ous parameters in equation (8) in a manner which 

poo 22 r 

identity \ e' c x dx = , one readily computes 

Ja -sc 


Consequently A = ^f g j + f g j j, the average of the initial and final velocities, and 
B = ^f g £ - When d n = 0 for all n the velocity is A at x = 0; and, inasmuch 

as the odd-order Hermite polynomials are odd functions of x, the presence of the terms 
containing these polynomials does not change the velocity at the origin. The even-order 


7 



polynomials, on the other hand, influence the velocity function in a symmetric (even) man- 
ner and, consequently, affect the velocity at the origin. 

2 2 

The nature of the exponential factor e _c x in the integral term of f g (x) and 

f (x) assures that f n (x) will be essentially flat outside of some neighborhood of x = 0. 

P r' 2 

Inasmuch as the terms of the summation each contain a factor of e -x , the neighborhood 

of x = 0 over which these terms alter fp(x) depends on the magnitude of c compared 

to 1. 

An example is shown in figure 1. Here an axial velocity distribution obtained by 
using only the two terms of fp(x) (eq. (3)) is compared with one obtained with the same 
values for the parameters except with d^ = 0.1. Thus, the initial and final velocities and 
the velocity at the origin are all unchanged, but the variation of velocity throughout the 
design region is radically changed. This distribution ^with dj = O.l) has appropriate 
characteristics for a contraction-cone design, that is, it is relatively short between the 
flat ends with smoothly varying curvature. Figure 2 shows a contraction cone designed 
from this velocity distribution together with several internal streamlines. Figure 3 shows 
the distributions of velocity along these streamlines. As expected the radial variation of 
velocity is greatest near the entrance where the curvature is greatest. 

A different kind of design velocity distribution is shown in figure 4. Here the initial 
and final velocities were prescribed to be 0.5 and 1.0, respectively, with the terms with 
Hermite polynomials all chosen to have zero coefficients except d Q = 0.6. Thus the max- 
imum velocity occurs at x = 0, where f g = A + dg = 1.35. Such a velocity distribution 
(one with the peak velocity between the ends) cannot be described with the original Tsien 
formulation. 

A duct designed from this velocity distribution is shown in figure 5 together with 
some streamlines, and the wall velocity distribution is shown in figure 4. The radial 
variation of velocity is noticeable at the minimum, where the curvature is relatively large. 
This result may be compared with that of figure 3 for the contraction cone where the rela- 
tively small curvature at the minimum results in a nearly uniform flow there. 

CONCLUDING REMARKS 

A method for generalizing the Tsien procedure of contraction- cone design has been 
presented. The class of design velocity distributions is enlarged in such a way that con- 
ditions far upstream and downstream are unchanged, and so that the derivatives required 
in the calculation can be obtained by a recurrence formula rather than by numerical dif - 
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ferentiation. The generalized method is no longer restricted to contraction cones but now 
permits the design of diffusers and converging-diverging ducts. 

Langley Research Center, 

National Aeronautics and Space Administration, 

Hampton, Va., September 26, 1972. 
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APPENDIX 


COMPUTER PROGRAM FOR LOW-SPEED DUCT DESIGN 

A computer program has been developed which will calculate the wall contour for a 
subsonic duct. The program is written in the FORTRAN IV language for use on CDC 6000 
series computers. Since it was desired to take an interactive approach to the problem, 
the program has been implemented on the LRC interactive graphics system using the 
CDC 250 Cathode Ray Tube (CRT). The program listing and a description of its input 
and output are presented in this appendix. 

Description of Program 

The program is basically divided into two parts. Part I of the program builds an 
1x4 design table where the stream function xj/, the axial coordinate x, the radial coordi- 
nate r, and the number of derivative terms N are column vectors; I is the number of 
row entries. The program computes the value of the stream function at any point x,r or 
the value of r for an arbitrary value of \f/ at some specified axial coordinate. By 
employing these two computations (each of which stores an entry in the design table) , the 
user can determine the neighborhood of the desired solution and approximate the bounda- 
ries of its convergence. 

Part II of the program computes the radial coordinates which agree with some spec- 
ified range of the axial coordinates and a fixed value of the stream function (streamlines) . 
In addition, it gives the corresponding velocity distribution in the duct and its axial and 
radial components. The streamlines are visually displayed on the CRT with a visual cue 
at the point where the velocity is no longer monotonically increasing. The plotting speci- 
fications are variable and may be input during program execution. 

Subprogram Index 

The following is an index of the subprograms called by this program and their 
sources. AUTHOR denotes routines written by the authors of this paper. CALCOMP 
indicates routines available as a part of the CalComp graphic output system. CRT indi- 
cates routines which are a part of the LRC interactive graphic system. LIBRARY denotes 
routines which are on the LRC computer complex system tape. The functions of the 
authors’ routines are also given. 
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APPENDIX - Continued 


Subprogram 

Source 

Function 

ABS 

LIBRARY 


AXES 

CALCOMP 


CALPLT 

CALCOMP 


CDC250 

CRT 


DAYTIM 

LIBRARY 


DERIV 

AUTHOR 

Computes the derivatives of f g (x), equation (6) 

ENCODE 

LIBRARY 


EXP 

LIBRARY 


FACT 

AUTHOR 

Computes the factorial of an integer 

FLOAT 

LIBRARY 


FOFR 

AUTHOR 

Evaluates ^ - f g (x) for routine ITR2 

FUNC 

AUTHOR 

Evaluates the integral in f g (x) 

IF IX 

LIBRARY 


ITR2 

LIBRARY 


LEROY 

CALCOMP 


LINPLT 

CALCOMP 


ME SAGE 

CRT 


MGAUSS 

LIBRARY 


NEXT 

CRT 




APPENDIX - Continued 


Subprogram 

Source 

Function 

NOTATE 

CALCOMP 


PARAMS 

CRT 


PNTPLT 

CALCOMP 


RECIN 

LIBRARY 


RECOUT 

LIBRARY 


RVALUE 

AUTHOR 

Computes r, velocities u and v, and G for ij/ at x 

SCREEN 

CRT 


SIGN 

LIBRARY 


SQRT 

LIBRARY 


STREAM 

AUTHOR 

Computes t// at any point x,r 



Program Input 

The first two cards should contain the velocity distribution function (free field). It 
will be printed as a part of the header on the first page of program output. 

The next input block should contain the velocity distribution function parameters and 
the parameters used in the iterative method to determine the radius of the duct. These 
variables should be input under the FORTRAN IV Namelist format. A description of these 
variables and the corresponding names used by the source program are as follows: 

Name 

Symbol 

Description 

$FPARAM 


Name required by input routine 

AG1 


Lower bound on the neighborhood of r 

AG2 


Upper bound on the neighborhood of r 

C2 

c2 

Velocity distribution function parameter 
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APPENDIX - Continued 


Name 


Symbol 


Description 


DELR 

D1 

D2 

EPS1 


EPS2 

MAXIT 

NACC 

VI 

V2 

$ 


^2 


N 


f g,i 


l S,i 


Initial size of the scanning interval r. ^ = r^ + DELR 
If r^ < r < r^ + p DELR = DELR/2 

Velocity distribution function parameter (see eq. (6)) 

Velocity distribution function parameter (see eq. (6)) 

Relative error criterion for determining convergence. If 


r i >s i> 


r i - r i- 1 


= 2^ implies convergence 


Absolute error criterion for determining convergence. If 
r^ = 2p | r^ - r^_ j j = E2 implies convergence 

Maximum number of iterations to be used 

Number of derivative terms to be used 

Desired initial velocity in the duct 

Desired final velocity in the duct 


Required by input routine 


The next input section forms the basis for the design table. Each card should con- 
tain an axial coordinate (columns 11-20, F10.4) and a radial coordinate (columns 21-30, 
F10.4). These coordinates should be chosen such that the stream function is specified 
throughout the entire field of interest. The value of the stream function is computed at 
these points and stored in the design table ordered on decreasing values of i p. 

The final input block is to be input at the CRT station. The variables in this block 
may be changed at any time during execution of the program affording interactive control 
over the program. By varying these parameters the user may take advantage of program 
options to (1) add and delete entries in the design table, (2) make limited changes to the 
velocity distribution function, and (3) vary the iteration scheme to achieve convergence. 
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APPENDIX - Continued 


Name 

Symbol 

Description 

AG1 


* 

AG2 


★ 

A1 

A 

Velocity distribution function parameter (see eq. (3)) 

A2 

2\j2cB 

Velocity distribution function parameter 

C2 


*• 

DELR 


* 

DPSI 


Increment from PSIMIN to PSIMAX 

DX 


Increment from XMIN to XMAX 

D1 


* 

D2 


* 

EPS1 


* 

EPS2 


* 

GDIST 


Length of Y-axis for velocity plot (in.) 

GDV 


Y-axis scale for velocity plot (units/in.) 

GMAX 


Maximum velocity computed 

GMIN 


Minimum velocity computed 

GOR 


Y-axis origin for velocity plot 

I 


Row number in design table 
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APPENDIX - Continued 


Name 

Symbol 

Description 

IP1 


Printing control for part I of program 

IP2 


Printing control for part II of program 

MAXIT 


* 

NACC 


* 

PSI 

i' 

* 

PSIMAX 


Maximum streamline to be computed 

PSIMIN 


Minimum streamline to be computed 

R 

r 

* 

RDIST 


Length of Y-axis for radial plot (in.) 

RDV 


Y-axis scale for radial plot (units/in.) 

RMAX 


Maximum radius computed 

RMIN 


Minimum radius computed 

ROR 


Y-axis origin for radial plot 

VI 


* 

V2 


* 

X 

X 

* 

XDIST 


Length of X-axis (in.) 

XDV 


X-axis scale (units/ in.) 

XMAX 


Maximum value of x for which \p is computed 



APPENDIX - Continued 


Name Symbol 


Description 


XMIN Minimum value of x for which \ p is computed 

XOR X-axis origin 

NOTE : 

The starred (*) variables are defined previously in the appendix. 

A sample input follows: 

F(X) = A1 + A2*I(SQRT(2PI)*E**(-C**2*X**2))DX + E ** - (X**2)*(D1*H0 + D2*H1) 

I = INTEGRAL FROM 0 TO X 

$FPARAM AG1 = 0.0, AG2 = 1.0, C2 = 1.0, DELR = 0.1, D1 = 0.0, D2 = 0.1, EPS1 = l.E - 6, 
EPS2 = l.E - 6, MAXIT = 200, NACC = 10, VI = 0.133, V2 = 1.0$ 


1.0 

1.1 

1.0 

.9 

0.0 

1.0 

-1.0 

.9 

-1.0 

0.0 

1.0 

.7 

1.0 

.6 

1.0 

.4 

-1.0 

.7 

-1.0 

.6 


Program Output 

On the first page of printed output, the velocity distribution function, its parameters, 
and the design table are printed. 

On the following pages, the radial distribution is shown for ip = PSIMIN to 
ip = PSEMAX incremented by DPSI over an axial range of X = XMIN to X = XMAX incre- 
mented by DX. The resultant velocity G and its axial u and radial v components 
are also included. 

The plotted output included the radial distribution curves, the velocity distribution 
curves, and the centerline velocity curve. They are displayed on the CRT during program 
execution with the capability of saving them for post-processing on the Calcomp plotter. 
The plot format is similar to that of the figures shown in the main body of the paper. 

Following is the printed output which corresponds to the input previously presented. 
Streamlines are shown for \p = 0.003, 0.006, 0.009, 0.012, and 0.013. The centerline 
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velocity distribution is shown at \p = 0.0. The proposed wall contour is streamline ip = 0.013. Finally, the source 
program listing is presented. 
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APPENDIX - Continued 

PROGRAM CMO.ONc < INPUT = 100 1 »0UTPUT = 1 00 L * TAPE5= I NPUT » T APE8=0UTPUT *TAP 
i ta ) 

LOW SPEED OUCT DESIGN 

DIMENSION PSI t ( l O'J ) » A<100>« Ra( 100)« MESG ( 1 0 ) » TONE C 8 ) ♦ XPLOT(500 
1)» YPl.OT (500) . NACUnfl) 

DIMENSION QPG r p ) « DM (?) * 01ST<?)» JRCD(2) 

DIMENSION 0ATc(2)t FcT<8«2) 

EUUI VALENCE (cOR»ORG ( 1 ) ) • (GOR,ORG (2) 1 i (RDV.OV(l))* <GDV«DV(?))» 

1 (WDIST»DIST <1 > 1 , (GDTST,DIST(2» 1 
COMMON /YQ]/ 0FLR»EPSHEPS2»MAX IT 
COMMON /YQ?/ MACCtPST «X,R,U« V«G 
COMMON /Y03/ A 1 ,A2«C?«D1 »Q2 
COMMON / Y 0 A / A G 1 • A G 2 

NAMEL I ST/FPAR am/ AG1, AO,2« C2* DEl.P* Dl» D2» EPS1» EPS2* 

1 M A X I T « MACC. VI ♦ VP 

PROGRAM INITIAL I /ATI ON 

CALL C 0 C 2 s 0 
CALL LEROY 

CALL SCREEN ( 1 . 0 ♦ l • 0 ♦ 0 . 8 > 

S2=SQPT (2.0) 

MN A= 1 0 0 
MPT=500 
1P1=IP2=0 
JNK 1 =0$ JNK2=- 1 
LPLT=0 

CALL D A Y T T M (oaTE» 

READ 1 1 0 » FCT 
110 FORMAT ( 8 A 1 0 ) 

READ (S'FPAPAM) 

C 

CALL “ARAMS 

CALL PARAMS Ol PShPSlI 

CALL PAPAMS ( 1 1.X * X ♦ 1( R »R ♦ 1LI » I > 

CALL “ARAMS (6LDELR,DELR,SLMAXtT»“AXIT) 

CALL PARAMS (ALtPSl ,FPS1 ♦alEPS?«E“S2) 

CALL PARAMS ( PL A 1 ♦ A 1 . ?L A? * A2 ) 

CALL PARAMS (aI.NACC»NACO 

CALL PARAMS ( PL AG 1 • AGl ♦ 3L AG? » AG? ) 

CALL PARAMS ( 3L I P 1 ’ I P 1 » 3L I P2 » I “2 ) 

CALL PARAMS < o | v 1 .VI .?LV?,V2,?LC2.C?> 

CALL “ARAMS ( PI..D 1 » O ] . ?LD? * D? ) 

CALL ME SAGE ( 1 . 3Rh TmSERT INITIAL VALUES FOR PSI EOIT»3b) 

CALL MESAGF (i.?Sm AMY FN KEY «FLE CON I I N' IF . 2S ) 

CALL NEXT (KEv) 

J3 = l 

GU TO bio 
120 CONTINUE- 
C 
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C PAPT I 

C GENERATE INITIAL PSI DI STRIBUT ION TABLE 

NA = 0 

130 Rt AD 149* X * R 

IE (EOF. 5) 180. ISO 
140 FORMAT ( 1 OX »2r 1 0 .4 ) 

1SI) NA = NA*1 

IF <NA.LE.MNA> GO TO 1 70 
PRINT 160. MNa.NA 

160 FOPMAT ( 3RH0M ax I MOM INITIAL PST DIMENSION EXCEEDED* 21 1 0 ) 

GO TO 130 

C COMPUTE PSI 

170 CALL STREAM 

PSIA(NA)=rSI'Sa(NA)=X‘»PA(NA)=R 
NAC (NA) = N A C C 
GO TO 130 

130 IF (NA.GT.MNA> NA=MNA 
J 1 =27 

ENCODE ( Jl * 19' .MESG) NA 
190 FURMAT (1X.I5.R1H INTTIAL PSI COMPUTED) 

CALL ME SAGE (l.MESG, I] ) 

C 

C 

C ORDER INITIAL PSI TABLE ( DESCENDING ) 

J3=NA-1 

DO 210 Jl=l*Ji 
J2=J 1 ♦ 1 S Jp= J l 
DU 200 J4 = J2 » niA 

IF (PSIM.I5) .( T.PSlA(J4n J5= J<+ 

200 CONTINUE 

IF (J5.FQ.J1) GO ID ?1U 

SAV = PSI A ( II ) <£PSI A ( Jl ) =P$I A < J5) SPST A ( JS) =SAV 
b A V = A ( J ) ) s a ( J l >=A( J5> S A ( J5 ) = S A V 
Sav = RA ( Jl ) 'fiRA ( Jl ) =RA ( JS) $RA ( J5» =SAV 
bAV = NAC ( Jl ) 1>M AC ( Jl ) =NAC ( J5) *NAC < J5) =SAV 

210 CONTINUE 
C 

C DISPLAY INITIAL PSI TABLE 

220 JRT=1 
230 J2=20 

J5 = 50 

DU ?6 0 J 1 = 1 • N a * J? 

J3=J1 * J2-1 
IF (J3.GT.NA) J3=DA 
DO 250 J4 = J 1 ♦ J3 

ENCODE ( J5.24-. .MES(i) J4 ,PSIA(Ja) .A(J4> .PA(J4> «NAC(J4) 

2a(j FORMAT (3p I = . ) 2 . 5H oSI = .Pl 3.6.3H X = ,FH.‘+.3 M R=«F8.4*3H N=.I2) 
CA|.L MFSAf.E ( l , Mf 5G. 15) 

2b j CONTINUE 

CALL MfSAGF ( i . 3 3 h AMY IN KEY I-ILL CONTINUE UISPLAY.33) 

CALL NEXT ( K F_ / ) 

265 CUNT INUF 
C 
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C EDIT INITIAL PS! TA6I E 

C 

CALL MFSAGE (1.25* END OF PSI RANGE DISPLAY, 25) 

GO TO (?7i).6?r'>* JPT 

270 CALL MFSAGE <1.24* RcGlN EDIT RANGE OP PSI. 24) 

CALL MFSAGE (l«?7H VARIABLES APE X, R. 1*27) 

280 CALL MFSAGE (1.37* EN KEY 1 WI|L COMPUTE P SI AT X AND R.37) 

CALL MFSAGE (1.38* EN KEY 2 WTlL INSERT PSI. X, AND R,36) 

CALL MESAGE (1.23H Fm KEY 3 WJlL DELETE 1*231 

CALL MFSAGE (1.39H EM KEY A Wll.L DISPLAY RSI, X, AND P AT 1.39) 

CALL MESAGE ( 1 • 32* FM KEY 5 WlIL DISPLAY PSI RANGE.*32) 

CALL MESAGE <].36H EM KEY 6 Wll.L END EDIT AND CONTINUE. 36) 

CALL MESAGE (1,37* Em KEY 7 Wil l. COMPUTE R AT PSI AND X»37) 

CALL MESAGE O.A]* Fm KEY 8 SImPL c TRANSFER TO PLOT R0UTINF*41> 
CALL MESAGE (1-38* FM KEY P WILL DELETE ENTIRE TA6LE.38) 

CALL ME SAGE <1.33* F,\i KEY 10 WILL COMPUTE Al» A2.33) 

CALL MESAGE (1.38* ANY OTHER F:J K* Y WILL DISPLAY OPTIONS, 36) 

C 

CALL NEXT (KEY) 

IF (KEY.EO. 1 ) GO TO ?90 

IF (KE.Y.EQ. 2) GO TO UG 

IE (KEY.EO. 3) GO TO 400 

IF (KEY.EO. A) GO 10 ^30 

IF (KEY.EQ.S) GO TO ??G 

IE (KEY.EO. 6) GO TO S4G 

IF {KEY.EO, 7) GO TO a4G 

if (kfy.eo.8) go ro g20 

IE (KEY.EO. 9) GO TO A9G 

IF (KEY.EO. 10) 00 TO GUO 
GU TO 280 
C 

C COMPUTE PSI 

290 CALL STREAM 

ENCODE ( 50 * 30c • MESG) PSI.x,R,NaCC 
300 FORMAT (5h PSt = .F 16.g,3h X=,E8.A»"<h R=,E8.4,3H N=,IA) 

CALL MFSAGE (1.MESG.S0) 

GO TO 280 
C 

C INSERT PSi IN lARLF 

310 IE (NA.FO.MNAi GO TO 390 

IF (PSI .LT.PST a ( 1) ) GO TO 320 
Jl = l 

GO TO 340 

320 00 330 J 1 = ? * N a 

IF ( (RSI .1 E.PSIA(J1-1 ) ) .ANO. (PSI.GF.RSJA ( Jl) ) ) GO TO IaO 
330 CONTINUE 
J I =N A ♦ 1 
GO TO 3*0 

340 00 35') J2 = J) , ' : « 
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J3=N A - J? ♦ J 1 

MSI A < J3*1)=PStA<J 3>*A (J3*l)=A( 13* t>RA ( J3* 1 > =9A ( J3> 

NAC ( J 3+ 1 )=NAC<J3) 

350 CONTINUE 

360 HSIA(.JI) =PSIS/> ( J1 ) = X6PA ( Jl) = R 
NAC ( Jl) =NACC 
N A = NA + 1 

370 ENCODE (76.38” .MF.SG) NA ♦ Jl . NAC ( J I > . PSI A ( Jl* ♦ A ( Jl ) »PA ( Jl) 

380 FORMAT (9h Total I=,T4.3H I = » 1 4 . 3h N=,I3«4X.5H PSI=.F16.6*4H X = .F 
1p.4,4H R= . F9 . 4 ) 

CALL MESAGE O.MESO.30) 

CALL MESAGE ( i .MESG(A) .47) 

00 TO ?80 

390 CALL MESAGE (a. 23" MAX I MAS RfFN REACHED. 23) 

00 TO 28 0 
C 

C DELFTF r»:T9r I IN PSI TAHLE 

400 UO 410 J 2 = I » N 

MSTA ( >2) = PSlA < 12+1 ) 

A( J2)=A( JP+1) 

K A ( J2 ) =R A ( J2+1 ) 

It A C ( J ? ) = N A C ( J ? * 1 * 

410 CONTINUE 
Nrt=NA-l 

ENCODE ( 30 . 4?i; • MESG) NA.I 

420 format (9w total i=.ts*I8.bh OFLETEO) 

CALL ME. SAG £ ( ’ . M E S G . 3 0 > 

00 TO 260 

c 

43C Jl=I 

GO TO 370 
C 

C COMPUTE » FOP PST AT X 

440 CALL pva'LuE ( t COOL * 

litCODF (50.45n.MESG) PSj.X 
450 FORMAT (5 h PS T = . f 1 6 . a . 5h X=»f16.8> 

CALL MESAGE (l.MESG.4?) 

ENCODE (5o.46o.MESG) R*G 
460 FORMAT ( 3H Ps.FIG.H.RH G=,Fl6.8) 

CALL. MESAGE (1.MESG.40) 
tNCODE (50.470.MESG) IJ.V 
470 FORMAT (3 h U= . f l 6 • b • 6H V=»fi6.h) 

CALL MESAGE ( 1 .MESG.41 ) 

ENCODE (50.480.MESG) NACC.ICODf 
480 FORMAT (3h H = . I ? . 1 2 M E*POP CODf=,I3) 

CALL MESAGE O.MESG.PO) 

CALL mfsAGF ( 1 . ?5h AMY FN KEY *i f L 1. CONTINUE. 25) 

CALL NEXT (KFy* 

GO TO 280 

490 NA = 0$PSl A ( 1 ) =r»SI A (2) =9.E10 

CALL. MESAGE ().14" TaRLE DFLETfD.14) 

Gu TO 2H 0 
C 

500 J3=2 

GO TO 5 1 0 
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c 

c 

r 

COMPUTE 

FUNCTION 

PARAMETERS 

Al 

42 

510 

C2 = ARS < C 2 > 






A 1 = 0. 5* <V1*V2) 

A2=S2*SORT (C2) « (V2-V) ) 

ENCODE ( 1 0 0 ’ 5P0 .MESG) VI . V2.C2. Al , A2 
52C FORMAT ( 1?H IMIT. VEL ,=»F8.4, 1?H FINAL VEL . = . F8 .4 » 4H C2= »F 1 0 . 6 , 6X . 
1 AH Al=,Fl6.8,4H Ak=,F16.8) 

CALL MESAGE <1,MESG,40> 

CALL MESAGE ( 1 .MESG (SM20) 

CALL MESAGE < I . MESG ( 7 ) ,40 > 

ENCODE (34. 53a, MESG) Dl.02 
530 FORMAT <4 h D1-.F13.6.4H D2=,F13.6) 

CALL MESAGE <l.MESG,34> 

GO TO (120.280) . J3 
C 
C 

C END INITIAL PS I EDIT 

540 Call MESAGE <3.?8h PfPEAT FN KfY a TO END EDIT. 28) 

CALL MESAGE ( ) . 38H AMY OTHER FN KFY WILL DISPLAY OPTIONS. 38) 

CALL NEXT (KEY) 

IF (KEY. NR. 6) GO TO P80 
IF ( IP1 .NF.0) GO TO HlO 
C 

C OUTPOT INITIAL PSI TABLE 

DO 60 0 J 1 = 1 . N a 

IF ( MOD ( J 1 — 1 . 35 ) . NE • 0 ) GO TO 505 
PRINT 550. n ATE ( 1 ) 

550 FORMAT ( 1H] //5(!X. "CONTRACTION CONE DESIGN TABLE*. 15X, "DATE*. 

1 5X. A i 0 > 

PRINT 560, r r.T 

560 FORMAT (/I OX.gOHVELOCITY DISTRIBUTION FUNCT I ON , 5X , 8 A 1 0/45X . 8A 1 0 ) 
PRINT 570, V), V2, Al, A2, C2 

570 FORMAT ( /?0X , j ! H IN I T . VEL . = , F9 . 4 ♦ 5X , 1 1 HE INAL VEL . = , F9 . 4 » 5X , 3H A 1 = ,F 
1S>.4,5X.3HA2 =»fo.4.5X,5HC**2=,Fs.2) 

PRINT 571.01, 02 

571 F ORM A T ( 53X , *0)=*, Fft.4, * D? = *. F8.4) 

PRINT 580 

580 FORMAT ( / /32X , 1HI , 17X.3HPSI , 19X » IhX, IPX. 1HP.5X, 1HN) 

595 CONTINUE 

PRINT 590. J1 ,PSI A ( Jn *A ( J1 > ,RA ( Jl> ,NAC( J1 > 

590 FORMAT ( 28X . I S . 3F20 .5 . 1 6 ) 

600 CONTINUE 
610 CONTINUE 
C 
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part II 

set-up plot parameters 

call PARAMS (4LXMlN,xMlN,2LDX,nX.4LXMAx,XMAX) 

CALL PARAMS (^LPSIMlN,RSIMIN*4LDPSI*nPSI»6LPSIMAX.PSIMAX) 

CALL PAPAMS ( OLROk ♦ RqR 4 3LRDV » Ro V » SLPDI ST* PD 1ST) 

Call PARAMS ( 3L GOR ♦ GOR ♦ 3LG0V » GOV » SLGD I ST ♦ GOI ST ) 
call PARAMS (Ol X0K,X0R«3LXDV»XnV.' : ;LXDIST»X0IST) 

CALL PARAvs („.lRMIN,PM1N.<*LKMAX,RmaX> 

CALL PARAMS (4LGMI'M»GMIN»4LGMAX,GmAX) 

XMIN = -4.(KDX = r. .4 SaMAx= 4.0 
RSlMlM=.0?7inPSI=i *0^PSIMAA=.0?7 
18=8 

JUATA=JPS=JGS=0 
CALL CALPLT (0. 0*0. 0.-3) 

TMA J=1 ,0$TMlN=0. I 
JXpCD=lHX*XHG T =4. J/1S.8SJNX=-1 

XUlST=l0.0iRDrST=i0.nsGOIST=10.0 
JBCD ( 1 > =1 hRSJpC0(2> =1 HG5>JN=*1 

jsym=osjS] ze = i < f.py=o. i 
mcurvf=i i 
c 

620 CALL MRS age ( 1 . 3 ?H Em KRY 1 W I ( L OISPLAY RSI RANGE , 32) 

CALL MESAGE (1.27H Em KEY 2 WI|_L COMPUTE DATA. 27) 

CALL MESAGE (l.46n VaPIARLES APE XMIN.DX, XMAX.PSIMIN.DPSI .PSIMAX.4 
lh) 

CALL MESAGE (1..37H EM KEY 3 WTlL SCALE PSI 0 1 STR I HUT ION, 37 ) 

CALL MESAGE <i,42h VaPIARLES APE XOR , XDV , XD I ST , ROR . PDV , ROI ST , 42 ) 
CALL MESAGE (I.37h Em key 4 WILL PLOT PSI DISTRIBUTION, 37) 

CALL MESAGE (1.35H EM KRY R WH.L SCALE G 0 1 STRlBUT ION , 35 ) 

CALL MESAGE (1.31H VaPIARLES ARE GOP, GDV . G0IST.31) 

CALL MESAGE (1.35H EM KEY 6 WT L L PLOT G D I STR I BUT I ON , 35 ) 

CALL ME STAGE < 1 , <*4M EM KEY 7 W 1 1 L DISPLAY NON-MONOTONE VELOCITY. 44) 

CALL MESAGE ( ) . 26H En KEY 8 WILL END PROGRAM, 26) 

CALL MESAGE <1.39* En KEY 9 k[i.L RETURN TO PSI RANGE EDIT, 39) 

CALL MESAgR <1,?Sh en KEY 10 WTLL NORM aL I ZE . 25 ) 

CALL MESAGE (1.38* AMY OTHER Em KEY WILL DISPLAY OPTIONS, 38) 

CALL NEXT (KEYl 


IF 

(KRY. ED. 1) 

GO 

TO 

G30 

IF 

(KEY. ED. 2) 

GO 

10 

640 

IF 

(KEY. F 0.3) 

GO 

TO 

mo 

IF 

(KRY. ED. 4) 

GO 

TO 

930 

IF 

(KFY.fO.5) 

GO 

TO 

1 090 

IF 

(KEY.EO.ft) 

GO 

TO 

930 

IE 

(KFY.E0.7) 

GO 

ru 

1 120 

IF 

(KRY.cO.R) 

GO 

TO 

1170 

IF 

(KRY.E0.9) 

GO 

TO 

970 

IF 

(KEY.EO.10) 

GO 

TO 

850 

uU 

TO 620 
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C 

C 

630 

C 

C 

640 


C 


650 

660 

C 


670 


680 


C 

C 


C 

C 


C 


690 

700 

C 


DISPLAY INITIAL PSI TA6LE 

JRT = ? 

GO TO ?30 

COMPUTE RADIAL DISTRIRUT ION 

RCURVE= (PSIMAv-PSIMIn) /DPS I 
JNK 1 = JNK 1 + 1 
LPLT=LPLT*1 

CHECK FOP MAXIMUM NUMBER OF CURVES 

IF (RCURVF.LE.pLOAT (MCURVE- 1 ) ) GO TO 660 
ENCODE (37.65-i.MESG) RCUPVt .MCURVr 
FORMAT (2X*Fl?.5»iIH MORE THAN ,Is.7H CURVES) 

CALL MES Ant U.MESG.37) 

GO TO 6?0 

RPT= ( XMAX-XMIM) /DX 

CHECK FOP MAXIMUM NUMBER OF POINTS 

IF (RPT.LF.FLPAT (MPT-l) ) GO TO 68.' 

ENCODE (37.67'..MESG) RPT.MPT 

FORMAT (2X.F1?. 5. llH MORE THAN .IS.7H POINTS) 

CALL MESAGE (A.MESG.87) 

GO TO 6?0 

ICURVE=IFTX (RrUPVt ♦] .5) 

NPT = TF I X ( RPT ♦ i .5) 

REWIND 18 S ,ITONE = 0 
RMAX=6MAX=-1 ,fq 
RM I f\|=GM T Nr 1 .Fo 
PS I =PS I M I N 


UO 790 J2= 1 » I rilRVE 
OT ONE = ~ 1 • F9 
IE0F=uSJ3=-l 
X = XMTM 
I CODE= 1 

DO 780 J 1 = 1 » Np T 

compute radius 

CALL RVALUE ( t CODE. ) 

IF (ICOOE.EO.n) Dp=R 

CALL PECOUT < ) A ♦ 1 ♦ IEOF * J? * J1 » PS I i* .R.U.V.G) 

IE (JTONE.NE.i-) GO TO 700 

CHECK EOp NON-MONOTONE VFLOCJTY 

IF (G.GT.GTONp) GO TO t>90 

TONE ( 1 > =J?iTOt.'P (?) =J1 6 TONE < 3 ) =pS I FT ONE (4 ) =X 

TONE (6) =Pf TONp (6) =U‘(>TONE <7) = V<bTONF (8) =0 

JTONE=l 

GO TO 700 

G T ONE = G 

CONT T NUF 

CHECK F Qp MINIMUM AND MAXIMUM YRLOT 


VALUES 


30 
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IF <R.LT.pMIN» RMIN=P 
IF (R.GT.PMAX* RMAX=P 
IF (G.LT.GMim GM1N=G 
IF CG.GT.GMAX> OMAX=g 
IF (IP2.MF.0) GO TO 7?0 
J3=J3*I 

C OUTPUT CURVE j? 

IF ( MOD ( J3 * 15 ) .NP. 0 ) GO TO 765 

PRINT 710. J?. D A 1 E ( 1 I » PSIM1N, PS1MAX, XM1H, XMAX 

710 format (ihi//,-»7x, iah r and g (tistrirution for curve .I3»20x.4hda 

1 Tt ,5 X. a 1 0//20 V .6HPSI = . F 1 4 . 6 . RX * ?Rj 0 * F 14 . 6 * 1 OX » 3HX= . F 14 .6 » 5X , 2HT 
20 , FI 4. 6 ) 

PRINT 7?G , RSr 

720 FORMAT ( SripPS T = * F 1 4 . 6 > 

PRINT 730 . V), V2 1 Al. A2» C2. MACC 

730 FORMAT (/5X, 1 l MINI T. VtL . = , F8 . 4 , 5X . 1 1 HF I NAL VEL .= ♦ F8. 4 , 5X * 3HA 1 = ,FR 
1 .4.5X. 3mAP=,Fo.4»5X,GhC**2=«F8.4»5X«?HN=» 13) 

PRINT 740, Dl • 02, LPLT 

740 FORMAT (4f)X,3or.l = ,F8.4*5X»3HD2s:,FM.4,SX*8HPL0T NO. ,15) 

PRINT 750 

750 F ORMAT ( / ! 3X, prPT , 14X , 1HX , 1 4. X , 1HR, 14 X » 1HU, 14X, 1HV* 14X, 1HG) 

765 CONTINUE 

PRINT 760, JI.X,R,U,V,6 
760 FORMAT <115*5<MS.b) 

770 COMTINUF 
C 

X=x+Dx 

IF (X.6T.XMAX) X=xM*x 
780 COMTINUF 
C 

PST =PS I +DR5 1 

IF (PSI .&T.PSt‘:aX) PSI=PSIMAX 
790 CONTINUE 
C 

C OISpLaY MINIMUM ANO MAXIMUM PLOT VALUES 

C 

t MCODE ( 80 « 80 . MESG ) RMIN.RMAX .GMIN.GMAX 
800 FURMAT (6 m Rm t M= , F 14 . 6 » 6H Rm a X= , F l 4 . 6 , 6H GM IN= , F 14 , 6 , 6H GMAX=,F14. 
16) 

CALL MtSAGE (1.16H MATA COmPijTEU.Io) 

CALL MESAGF (1.MESG.4Q) 

call mesage ( r » mes6(5) «40» 

O O A T A = ) 

GU TO 6?0 
C 
C 

6 10 IF ( JO A TA.NE.ir> GO TO 820 

CALL MESAGF <7.29HCAMMOT SCALE NO DATA C0MPUTED*29) 

GO TO 6?0 
C 

C COMMUTE x AM) R SCAI..FS 
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620 AOR=XMIn 

A D V = ( X M A X - XM I m > / X i) I ST 
HORsRMIN 

ROV = < PMAX-RM In) /ROI ST 
ENCODE ( 3ft .83" .MESG) XOR.XDV 
630 FORMAT (5 h XQo= , F 1 2 . S , 7H X0V=,Fl2.5) 

tNCODE <3ft.84'‘,MESG(ft> >R()R*RD\/ 

840 FORMAT (5H R09=,F12.5,7H RDV=,E)2.S) 

CALL MESAGE (t, MESG. 3ft> 

CALL MESAGE (l .MESG (ft) *36) 

JRS=1 
GO TO 620 
C 

c NORMALIZE last monotonf stream LINE 

850 JPLT=1 

IF (JNKl.NE.JfK2) GO TO 920 
UO 900 J 1 = 1 • NPT 

ARLOT ( Jl ) = X P L G T ( J 1 ) /YPLOT (NPT) 

YPLOT ( J1 ) r YPLOT < Jl> /vPLOT (NPT) 

IF (MOD ( Jl-1 05) .NE.O ) GO TO ««0 
PRINT HftO , OATF ( 1 ) 

860 FORMAT ( 1 H 1 // ?7X , 46HN0RM AL I ZE D CURVE FOR LAST MONOTONF STREAM LINE 
I ♦ 1 OX » AHDATE «5x . A1 0 > 

PRINT 720. PSI 
PRINT 870 

870 FORMAT (//33X.2HPI » 19X* IPX* 19X, 1HR) 

880 PRINT 890. Jl . x.PLOl ( II )♦ YPLOT ( Jl ) 

890 FORMAT ( 2g X . I ) o » 2F20 . 5 ) 

900 CONTINUE 
C 

tNCODE (80.91(.. MF.S6 ) X.PLOT ( 1 ) .XPLOT(NPl) .YPLOT (NPT ) .YPLOT ( 1 ) 

910 FORMAT (<*H X1-.F13.6.4H X 2= , F 1 3 . 6 • 6X . AH Rl=»F13.ft«4H R2=»F13.ft) 
CALL MESAGE (). MESG. AO) 

CALL MESAGE ( l .MESG (S) «3A) 

CALL MESAGE (1.36H INPUT SCALE FACTORS PRESS ANY KF Y ,36) 

CALL NEXT (KEY)) 

JNK2= JNKl 
920 CONTINUE 
GO TO 950 

C READ p LOT DATA JPLT=1, RPLOT, JPL0T=2.GRL0T 

C 

930 JPLT = KE Y/S+ 1 

IF ( ( JPLT.FO.) ) .AND. (JRS.EO.O) ) GO TO 940 
IF ( ( JRLT.EU.p) . AND. ( JGS.EO.O) ) GO TO 940 
GO TO 950 

940 CALL MESAGE (l.39*-i CANNOT PLOT. DATA HAS NOT BEEN SCALFJD.39) 

GO TO 620 
C 

950 CALL AXES (0... .0.0,0. O.XOIST .XOR, XDV.TMA J.TMJN, JXBCO. XhGT.JNX) 

call axes ( 0 . m . 0 . vJ ♦ 9 a . 0 , D I ST (J»LT> .ORG(JPLT) .DV(JPLT) . TM a J. TM I N , JR 
ICO ( JPLl ) . XHGT. JN) 
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REWIND IB 

UO 99u J2 = l»IriiRVK 

IF ( KEY .£(•>. '0 ) GO TO 970 

UO 960 11 = 1 »NPT 

call pec in (i^.i.ic.ncv.jpt.pst.x.r.u.v.g) 

aFLOT ( J1 ) =X 
YPL0T(J1>=9 

I F (JPLT.FQ.?) YPLOT I Jl ) =G 
9b C CONTINUE 

C 

970 aFLOT (MuT + 1 ) =x()BY>aPLOT (NPT + ?) =xOV 

YPLOT (NPT+l )=r>PG( JPLT) * YPLOT (NPT.R) = DV < JPLT ) 

C PLOT DATA 

call limpet (xrlot .yplot.nrt. i .jsym.jp. jsize.oj 

c 

C NOIaTf PLOT 

JJ=17 

encode < J3.9rk.mesg> psi 

980 FORMAT (4MPS1=,F13«6) 

HHT = 0. 1 

PY = FLOAT ( J2-1 ) * (PhT* ) .2) 

RA1=?.S 

Call notate <-pxi »py.pht.mesg»o.o« j3> 

PA = (b.0/7.0) “FI.OA I ( JO) *PHT + 0. 14-PX1 
RY=PY+0.0 a 

CALL PNTPLT (pX.PY.J?,!) 

IF ( *EY .EU. 10 > GO TO 991 
990 CONTI NUF 
991 CUNTINUF 

encode <is. i.ooo»MbSG> vi 

1000 FORMAT ( 3m V 1 = . F 1 3. 6) 

HY = DI ST ( JPLTl -° hT 

Call notate <-PXl .PY,PHT.M£SG.0.0. 16) 

ENCOOF ('U.10mT*MESG> VP 
loio format (3hv2=. fi3.6> 

PY=PY-2.0^PPT 

call NOTATE (-PX1.PY.PhT.MESG. 0.0. 16) 

encode ( i*. io?g*mesg> c? 

1020 FuRMAT <5hC**?=.F13. a) 

PY=PY-?.0opHT 

CALL NOTATE (-PX1.PY.PHT.MESG. o.o. 18) 

ENCODE (1a.10to.MESG) Dl 
1030 FORMAT (3 md1=.F1 3.6) 

RY=PY-2.0»PHT 

CALL NOTATE (-PX1 .PY.PHT.MfcSG.O.0. 16) 

ENCOOE < lS.lO'*0»MtSG> U2 
I jAO FORMAT ( 3 ho2= . f 1 3. 6) 

H y = py-? . 0*ohT 

CALL ‘-IOTATE (-PX1 .RY.PHT.MESG.t.O. 16) 

ENCODE ( Is. 10 f-).MF SGi A1 
1 C 5 0 FORMAT ( 3mA1 = .f1 3.6) 
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HY=PY-2.0*PHT 

Call notate (_PX1.PY,PHT»meSG.0.0,16> 

ENCODE ( 1 (S * IOaO.MESG) A? 

1060 EOomat OMA?=.ri3.6i 
PY=PY-?. O^PHT 

CALL NOTATE (-PX1 .PY.PHT.MESG.0,0. 16) 

ENCODE (1A.1070.MESG) LPLT 
1070 FORMAT (9MPL0T NO. ,TS) 

PY=PY-2.0*PHT 

CALL NOTATE (- U X1 .PY.PHT.MfSG.n.O. 14) 

C 

IF (JTONE.EO.o) GO To 10RQ 
IF (KEY.EO.10) GO TO 1080 

C FLAG NoN-MONOTONE VELOCITY ON PLOT 

PX= (TONE (4) -XnP) /XDV 

PY= (TONE (3*JPl T+2) -OPG ( JPLT) ) /nV(.lPLT) 

CALL PNTPLT (ox,PY.n,3) 

1080 CALL CAlPLT ( l 4 . 0 » 0 . n » -3) 

GO TO 6?0 

1090 IF (JOATA.NE.o) GO To 1100 

CALL MESAGE (7.P9HCANN0T SCALE.NO DATA COMPUTED*?'-') 

GO TO 620 

C COMPUTE G SCALES 

1100 GOR=GM I N 

GD V = ( GM.A X -GM I m > /GD I S T 
ENCODE (3a.Hic.MESG) GOR.GDV 
1110 FORMAT (5h G0p=.F12.S.7h GDV=.F12.5) 

CALL MESAGE (I.MESG.GS) 

JGS=1 
GO TO 620 

1120 IF (JTONE.NE.r-,) GO TO 1130 

CALL MESAGE (1.23m VELOCITY IS MONOTONE. 23) 

GO TO 620 

C DISPLAY NON-MONOTONE VELOCITY 

1130 ENCODE <43*llA0.MESG) TONE < 1 > . TONf < ?) 

1140 FORMAT (7h CUp\/E *F3.Q«5H PT.E4.0.22H NOM-MONOTONE VELOCITY) 

CALL MESAGE (Y.MESG.aI) 

ENCODE (50.llco.MESG) TONE (3) .TONL (4) , TONE (5) 

1150 FORMAT ( 5H PSt=.F12.S.5h X=.E12,S.4H H=.F12.5) 

CALL MESAGE (1.MESG.C0) 

ENCODE (50 . 1 lAO.MESG) TONE (6 ) ♦ TONE < 7 ) . TONE ( H ) 

1160 FORMAT ( 3h U=.ri?.5,SH V=»Fl?.S.6H G=.F12.5) 

CALL MESAOE (I.MESG.eO) 

GO TO 620 

NORMAL PROGRAM STOP 
1170 CALL MESAGE (a.31H REPEAT FN KpY A TO END PROGRAM. 31) 

call MESAGE (A.JftM AMY OTHt R F-l Key WILL DISPLAY OPTIONS, 38) 
CALL NEXT (KEY) 

IE (KEY.NE.8) GO TO a?0 
RR 1 NT 1180 
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11P0 FORMAT <2SHl NORMAL FOD OF JOB ) 

STOP 

end 

BURPOUTiNF stream 

COMPUTE VALUE OF The STREAM function psi at X AND « 
X = AXIAL COORDINATE 

R = RADIAL COORDINATE 

COMMON /TOP/ MACC»F>ST «X,R 

COF( J)=(-1 .0)«*( J-l)*R**(2*J)/(2.o**JP*FLOAT( J)*F4CT( J-l)*»2) 

BSI =0.0 

00 10 Jl=) • NArC 
Jc! = 2*.II-2 

PS T =P$ I ♦ COF ( J I > *nt.RTV ( J2> 

CONTINUE 
RSI=0.S*PSI 

RETURN 
END 

surrouttnf RVALUE (I CODE l 


G 

U 

V 


COMPUTE «♦ U. Vi AND r, FOP 
IS THE RESULTANT VELOCITY 

IS THE AXIAL component 

IS thf radial component 

USE INTpRVAL-HALVING MPTHOD 


PSI 


AT 


X 


COMMON / Y D 1 / nFLRiEPSl *EPS2,MAYIT 
COMMON /YD?/ NACCiPSI *X,R,U»V.r, 

COMMON /YOA/ AG),AG2 
EXTERNAL' fofr 

COF ( J) = (-1 .0) J-l) »P**J2/ (2.0*° <?*J-3> *FACT ( J-l>**2> 
COFF(J)=(-1.0 i**(J)*p**J3/(2.C^*(P*J-?)*FLOAT<J)*FACT(J-1)**2) 


CALL I T R2 <R.aG 1»AG2,UELP.F0FR,FPS1,FPS2.MAXIT,IC0DF> 
IF < ICOoE.FU.'M GO TO lo 
K = IJ=V=G=-Q09.'i 
RETURN 


U=V=0.0 

DO AO Jl=liNArC 

J^=2*JI-? 

J J=2* J 1 - I 


IF ( ( I2.NS .0) .DR. ( AflS (R) .GI . 1 .P-IP) ) GO TO 20 
U=2. O^DFRl V ( JP> 

GO TO 30 
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CONTINUE 

U=U*C0F<J1)*DfPIV(J?> 

CONTINUE 

V = V + COFF ( JJ ) *0,FRI V ( J3) 

CONTINUE 

U=0.5*U 

V=0.5»V 


G=SO«T <U**?*V**2> 

RETURN 

END 

function FACT <J> 

factopjal function 
fact=i .0 

IF (J.EO.O) Re TURN 
IF (J.F0.1) RfTIJRN 

DO 10 Jl=?«J 
FACT=FACT*ELOaT ( J l > 

CONTINUE 

RETURN 

END 

FUNCTION FOFR <R) 

Evaluate function fur it--? 
COMMON /Yn?/ m a CC « PS t « X * PR 


S AV 1 =PS I 
KRsR 

CALL STREAM 

fofr=r$i 

PSI=SAV] 

FOFR=PS i -fqfr 
RETURN 

End 

function oeRiv <J) 

COMMON /YD?/ MACC»PSI*X 
COMMON /Yn3/ A I ♦A?»CR«D1 

DIMENSION HERm(50>» hF w ( SO ) » F OF X ( 1 ) » ANSU) 
EXTERNAL func 

COMPUTE JTH DERIVATIVE OF VFLOC1TY 
SR=SORT (2.0*3. 1 4lSv?OS> 

EX=E'XP(-Cp*X*o?)/SR 

E X 1 =E <P < -X i 

HE PM U ) =E x 

HE PM (?) =-?.0«r ?*X*HFRM ( 1 ) 

HEP(l)=1.0 

HER (?) =?.(I*X 

HER ( 3) =4. 0 


FUNCTION 


APPENDIX - Concluded 


HEP <4) =R.O"X"t -)-12.0"X 
C 
C 

IF (J.NF.O) Go TO 20 
PL1=0.0T ; RL2 = X 
IF (PL?.GT.PLi ) G'l TO 10 
PL 1 =XSkL2 = 0 . 0 
10 CONTINUE 

CALL MO A US S <pL1»RL?.}0»ANS»FUuC»FOFX.1) 

DEPI V = A1 ♦STGN n ,0»X> *A2*ANS*EX 1* (nl"HEP ( ) ) ♦ I)2 "hER (?) ) 
PL TURN 
C 

20 IF ( J.NF. ] ) G') TO 30 

OEPI V=A?"FX-E> 1 " ( I > 1 "mF K ( ? ) ♦D2"hER(3> > 

RETURN 

C 

30 IF (J.NP.2) GO TO 40 

U£PI V = -?. 0*A2"C ?"a"EX*EX1« <(.U"HFR (3) ♦^"HEP (4) ) 

PE TURN 
C 

A0 her ( J* n =?.0*y*HF" < J) - 2 . 0 "FLOA r ( j-l ) *Hb'P ( -J— 1 ) 

HEP < J + 2) = ?.0*x"HF.m < J*1 > -2.0 "FLOAT ( J) "HER ( J) 

C 

iF (J.FO.O) Go TO 50 

HE PM ( J-l ) =-?. o"C2" ( X * H E R M ( J — 2 ) ♦ FLO AT ( J-3> "HFPM ( J-3) ) 
50 HEPM ( J) =-?.0*o?* ( A»nrRM ( J-l ) +F| OAT ( J-2 ) "HFPM ( J-2 ) ) 

51=1.0 

IF (MOi.)( J.P) .-if .0) SI =-1.0 

UEPI 0 = A?"4FPM f J) *Sl*FXl* (01 "HEP < J* 1 > ♦0?*HFR ( J + ?> ) 

RETURN 

END 

SUBROUTINE FUnr (X.FOFX) 

UIMFNSION FOFx ( 1 ) 

COMMON /Yn.V Al.A2.CP 

C EVALUATE VELOCITY FUNCTION FQP MGAIJSS 

bR = SOP T (2.0*3. 1 4150255) 

FOFX=FXP ( -C2" y ""2 ) /SP 

KETUPM 

END 
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Figure 1.- Axial velocity used for contraction-cone design ^dj = O.l) compared with that 
obtained with di = 0. Nonzero parameters: A = 0.5665; B = 0.4335; c = 1. 




Axial coordinate, x 


Figure 2.- Coordinates of wall contour and some streamlines 
for design velocity of figure 1. 
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Axial coordinate, x 

Figure 4.- Design (centerline) and wall velocity distribution for duct with an area 
minimum. Nonzero parameters: A = 0.75; B = 0.25; d n = 0.6; c2 = 0.5. 




Axial coordinate, x 

Figure 5.- Wall contour and some streamlines for the design velocity of figure 4 
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